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Using a matched asymptotic expansion we analyze the two-dimensional, near-critical 

'"JT ' reflection of a weakly nonlinear, internal gravity wave from a sloping boundary in a 

fT , ! uniformly stratified fluid. Taking a distinguished limit in which the amplitude of the 

incident wave, the dissipation, and the departure from criticality are all small, we obtain 

a reduced description of the dynamics. This simplification shows how either dissipation 

or transience heals the singularity which is presented by the solution of Phillips {The 

Dynamics of the Upper Ocean, 1966) in the precisely critical case. In the inviscid critical 

^ ' case, an explicit solution of the initial value problem shows that the buoyancy pertur- 

^Tj I bation and the along-slope velocity both grow linearly with time, while the scale of the 

reflected disturbance is reduced as 1/t. During the course of this scale reduction, the 

stratification is 'overturned' and the Miles-Howard condition for stratified shear flow 

(^ • stability is violated. However, for all slope angles, the 'overturning' occurs before the 

0^ [ Miles-Howard stability condition is violated and so we argue that the first instability is 

convective. 

_^ , Solutions of the simplified dynamics resemble certain experimental visualizations of the 

jy^ ' reflection process. In particular, the buoyancy field computed from the analytic solution 



is in good agreement with visualizations reported by Thorpe & Haines (1987) J. Fluid 
. ■ Mech. 178, 299-302. 

Oj! One curious aspect of the weakly nonlinear theory is that the final reduced description 

is a linear equation (at the solvability order in the expansion all of the apparently resonant 
nonlinear contributions cancel amongst themselves). However the reconstructed fields do 
contain nonlinearly driven second harmonics which are responsible for an important 



H I symmetry breaking in which alternate vortices differ in strength and size from their 

immediate neighbours. 



1. Introduction 

1.1. Reflection of internal waves on a sloping bottom 

Understanding the intensity and spatial distribution of turbulent vertical mixing in the 
ocean is an important problem in physical oceanography. Ocean models require accurate 
parametrizations of turbulent mixing to make realistic predictions of the transport of 
heat, salt and chemical species. Because the ocean is stably stratified, vertical mixing is 
inhibited and convection to great depth occurs only in restricted high latitude regions. 
But the fluid which has reached the abyss by convection must ultimately return to the 
sea surface so as to maintain a quasi-steady state. How and where this return flow occurs 
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remains obscure to this day. Munk (1966) and Armi (1978) suggested that significant 



vertical mixing takes place at sloping boundaries. Indeed, a recent report by Polzin 



et al. (1997) shows greatly elevated mixing rates above rough topography, especially in 



the deepest 150 m of the Brazil basin. 

bandstrom (1966J was the first to propose the oceanic internal wave field as a possible 
source of the energy which is needed to activate strong mixing near sloping boundaries. 
Internal waves are well documented in nature (see Munk 1981) and they obey simple, 
but unusual, reflection laws at a rigid boundary. In optics or acoustics, the incident 
and reflected wave rays make the same angle with respect to the normal to the reflecting 
surface, whereas for internal waves they make the same angle with respect to the direction 
of gravity. When internal waves reflect off a sloping bottom, the reflected wave has the 
same frequency to as the incident wave. But because the wave frequency is related to 
the direction of propagation by to — N sin f3 preservation of lu implies preservation of the 
angle /3 (N is the buoyancy frequency and, as shown in figure |l|, (3 the angle between 
the group velocity and the horizontal). This reflection law leads to a concentration of 
the energy density into a more narrow ray tube upon reflection. A fraction of the energy 
of the incident wave is radiated away as the reflected ray, but a part of the energy is 
lost due to turbulent dissipation, the rest being converted into irreversible mixing which 
accounts for the vertical transport of mass and matter. 

Probably the most effective situation for boundary mixing arises when an oncoming 
wave reflects from a bottom slope which nearly matches the angle of wave propagation. 
At this critical angle, the analytic theory of internal waves reflecting from a uniformly 
sloping bottom (Phillips 1966) predicts that the reflected wave has infinite amplitude and 
infinitesimal wavelength. These unphysical results signal the failure of the idealizations 
(for example, linear waves and inviscid fluid) made by Phillips (1966)|. 



bandstrom (1966J reported observational evidence of energy enhancement at a par- 
ticular frequency. Later Eriksen (1982 and 1985) presented data showing energy and 
shear enhancement near the critical frequency at a few mooring sites. Eriksen inter- 
preted deviations from a linear and inviscid theory as evidence for dissipation through 
shear instability and nonlinear interaction. The inviscid nonlinear case was theoretically 



considered by Thorpe (1987) who showed that singularities may occur at other frequen- 



cies when a finite amplitude incident wave interacts resonantly with its own phase locked 



reflected wave. [Gilbert (1993)| has also studied this phenomenon on the continental rise 
and slope off Nova Scotia. Energy enhancement at the critical frequency was significant 
at the 95% level for 8 of the 30 tests performed even though the overall concavity of 



the slope should have slightly inhibited the phenomenon. Finally, Eriksen (1998) has 
recently reported striking observations made on the steep flank of a tall North Pacific 
Ocean seamount. Eriksen found a dramatic departures from the quasi-universal Garett- 
Munk spectral model near the bottom in a frequency band centered on the frequency for 
which ray and bottom slopes match. 

1.2. Laboratory and numerical experiments 



Following the first observational studies, Sandstrom (1966) performed a laboratory ex- 



periment that clearly demonstrated the amplification that results from internal wave 



reflection off a sloping bottom. Then, Cacchione & Wunsch (1974) showed that at the 



critical angle, the strong shearing motion becomes unstable and results in the formation 
of a serie of periodic vortices. Overturning of these vortices produces mixed fluid that 
propagates into the interior as regularly spaced layers all along the slope. 

i'horpe & Hames (1987J reported evidence of the three dimensional structure of the 
boundary layer but they were unable to reproduce the formation of the vortex array 
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Table 1 . A summary of experimental parameters. 7 is the angle between the bottom slope and 
the horizontal. Nsmff the frequency in s~^; h is the depth of water in centimeters, ao is the 
ampUtude of the motion of the paddle in centimeters; Re=(^^ N sin /3 / u is the Reynolds number 
defined by Thorpe & Haines (1987) and used by Ivey & Nokes (1989). Re = N sin 13/ i^K^ is 
the Reynolds number defined by Slinn & Riley (1996) and us. a is the nondimensional param- 
eter defined in (2.12). e = (atana)^'^ is the small parameter used in our expansion and g is 
dissipative parameter defined in (3.8). 



seen by pacchione fc Wunsch (1974)|. The absence of vortices were also noted by Ivey 



& Nokes (1989) when they studied the mixing efficiency in the case of the breaking 



of monochromatic incident wavefield uniformly distributed over the sloping bed. Taylor 



(1993) was particularly interested in the deca ying turbulence and the mixing generated 
in the boundary layer. De Silva et al. (1997) considered the case of a small ray's width 
compared to the bed length. 



Ivey & Nokes (1989), |De Silva et al. (1997)| and [Mac Phec (1998)| have also shown 



the formation of filaments by instabilities as waves approach the critical frequency by 
the action of the incident and reflected waves. Exactly at the critical angle, instead 
of producing fine structured filaments, the waves produced turbulent vortices at the 
boundary. These vortices apparently mix boundary fluid which is presumably expelled 
along the isopycnal corresponding to the new density of the mixed fluid. As suggested 
by [Caldwell et al. (1978) , these intrusive layers could also explained the presence of a 
highly "stepped" temperature profile as the steep slope is approached at Lake Tahoe, 
California. 



On the numerical side, Slinn & Riley (1996) have shown the creation of a thermal 
front moving upslope at the phase speed of the oncoming wave in the turbulent boundary 
layer. For a steep front, the thermal front resembles a turbulent bore exhibiting nearly 
continuous localized mixing, whereas for shallower slope, the mixing is observed across the 
breadth of the domain and is temporally periodic. The internal wave field continuously 
pumps fresh stratified fluid into the mixed layer, while simultaneously extracting the 
mixed boundary fluid. 31inn &: Riley (1998b) have pursued their published work and have 
recently reported very detailed results on the effects of the slope angle on the turbulent 
layer thickness, mixing period and mixing efficiency. An advantage of the numerical 
approach is that it is possible to perform experiments over the shallow bottom slopes 
that are more typical of oceanic conditions. 



Other recent numerical work of Javam et al. (1997) showed that for near critical waves, 
the instabilities were triggered at the bed, while for supercritical waves the instabilities 
develop away from the bed. Javam et al. (1997)] also showed the nonlinear creation of 
harmonics. 



4 T. Dauxois and W.R. Young 

1.3. Questions and strategy 
Our goals in this paper are to understand the role of the nonlinearity in the reflection 
process, and to to characterize the instabilities responsible for the transition to turbulence 
near the slope. We also emphasize the role of transience and dissipation in healing the 
singularity which occurs when the reflection is critical. 

The paper is organized into four different sections. The formulation of the problem 
is described first in Section 0. We also present the main effect of the reflection process 
in the outer region. In Section 0, we focus our study to the reflection problem in the 
boundary layer close to the slope. In Section 0, we derive the explicit solution for the 
inviscid case, whereas the viscous effects are presented in Section 0. Finally, Section pi 
contains the summary and conclusions. 



2. Analysis 

2.1. Formulation in the coordinate system of the slope 

We will consider a two-dimensional, nonrotating, incompressible Boussinesq fluid, with 
constant Brunt- Vaisala frequency N and a uniformly sloping bed (angle of inclination is 
7) shown in figure nl We do not make the approximation of small inclination. In lakes, 
bed slopes are in the range 2°-20°, whereas in the ocean the r.m.s. slope of the sea bed 
is roughly 4° . In the vicinity of seamounts slopes can be considerably higher in the range 



6°-26° (see pe Silva et al. (1997)| and references therein). And most experiments use even 
larger angles of inclination. 

The incident wave is a nearly monochromatic group of internal gravity waves. The 
angle between the incident group velocity and the horizontal is /? and thus the carrier 
frequency of the group is A^ sin/?. 

Our analysis employs a slope-oriented coordinate system in which x is the distance 
along the slope and z is the distance normal to the slope. In terms of these tilted coor- 
dinates the stratification of the Boussinesq fluid is 

p = Pq[1- g^^N'^{zcos^ + xsin^) - g^^b] , (2.1) 

where b{x, z, t) is the buoyancy perturbation of the resting stratiflcation. 

Denoting by (u, w) the components of the velocity field, the equations of motion are 

Du 

siy\^b + px ^ vV'^u, (2.2a) 

— C0S7 6 + P2 = i^V^w, (2.2b) 

Db „ „ „ , ^ 

— + A^sin7u-|-iV^cos7u; = kV^6, (2.2c) 

Ux + w,^Q, (2.2d) 

where D/Dt — dt + udx + wd^ is the convective derivative. 
The boundary conditions at z = are 

u^w^O, N^cos-f + b,, = 0. {2.3a -c) 

The condition in (|2.3|c) is that there is no diffusive flux of buoyancy through the slope. 



Wunsch (1970) and Phillips (1970) have shown that this condition at a sloping boundary 
produces a flow near the wall in a stably stratified fluid. Under laboratory conditions 
using dissolved salt this flow is very small (but detectable). This no flux condition will 
cause some complications later. 
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2.2. The incident wave and nondimensional variables 

The situation we envisage is shown in figure 0. The incident wave nearly satisfies the 
critical condition that 7 = /3. Consequently the group velocity of the reflected wave 
makes a very shallow angle with the slope. The figure also anticipates some aspects of 
our analysis: when the reflection is nearly critical the advective term becomes important 
in a region near the slope. This quadratic nonlinearity results in the production of a 
second harmonic. This nonlinearly generated wave can then radiate away from the slope 
(see Thorpe 1987). The angle 6 with which the second harmonic leaves the slope is not 
necessarily shallow. 
Using the convention 

^=-■02, w = 'ipx, (2.4) 

for the streamfunction, we consider an 'incident wave train' 

[ip.b] sa [1, NK]A{t - sz) cos{K sin a x + Kcoaa z- Nainpt) (2.5) 

impinging on the slope. We have introduced the wave number K and, as suggested by 
figure 0, 

a = /3 + 7. (2.6) 

The envelope, A in (2_^), changes slowly relative to the space and time scales of the 
carrier wave; s~^ is the vertical group velocity of this slowly modulated incident wave. 
Thus, near the slope, the incident wave switches on slowly with a prescribed amplitude 

m- 

W e now nondimensionalize (2^) using the scales suggested by the incident wave in 
( ^.5[ ). The envelope can be written as 

A = 1pn,axA, (2.7) 

where tpraax is the maximum amplitude of the streamfunction and so max(A) = 1. The 
nondimensional variables are 

[x,z] = K[smax,cosaz], t = NsmPt, (2.8a-c) 

['0,w,w,6,p] = iprnax['>P, K COS a u, Ksuia w, KNb, Np\. {2.8d-h) 

Using the variables above, the nondimensional equations of motion are 

Du tana sin 7 , 1 „, 

Dt smp smpcosa Re 

Dw cot a cos 7 , 1 „9 ,„ „, , 

-jT- + —-.Pz - . ^ : b = ^V^u;, (2.9b) 

Dt smp smpsma ne 

Db sin 7 cos a cos 7 sin a 1 „n, ,„ „ , 

-— + '^ — u+ '^ — w = — V^fe, 2.9c) 

Dt sin/3 sin/3 Pe ^ ^ 

u^ + w^= 0. (2.9d) 



The differential operators in (2.9) are 



V^ = sin^adl+cos'^adf, (2.10a) 

^^=dt + aiud., + wd,). (2.10b) 



In (2.9), we have dropped the "hats" which would otherwise decorate the nondimensional 
variables. The dimcnsionlcss dissipation parameters are the Reynolds and the Peclet 
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numbers defined as follows: 



TV sin/3 , „ TV sin/5 ,„ , , ,, 

i?e= — —f- and Pe ^ — —f-. 2.11a,6 



The other nondimcnsional parameter in (2.1(]|b) is a measure of the nonlinearity 



... . . V'ma.. 2.12 

2N sm p 
Thorpe (1987)| gives a useful physical interpretation of the amplitude parameter a: the 



maximum slope of the isopycnals in the incident wave (2.9) is atan/3/(l — a) if a < 1 
or infinite if a > 1. In the sequel we will deal exclusively with the weakly nonlinear 
case in which a is small. The parameter a is also related to the internal Froude number 
fr= max(u2/A'^) by 

sin 2/3 cos ^ 

Fr=a ^— . (2.13) 

sm za 

Some authors prefer to use the minimal Richardson number Ri = Ff^ . Thus, with the 

appropriate geometric factors, a is simply related to the other measures of nonlinearity 

used in the literature. 

In this nondimcnsional and slope-oriented coordinate system, the dispersion relation 

of an inviscid linear internal wave (with all fields proportional to exp[i(fca; + mz — tot)]) is 

k sin a cos 7 — 771 cos a sin 7 
0J^± \ J=, (2.14) 

sm p^y [K sm a)^ + [m cos a)^ 

with the corresponding group velocity 

/ du) du) \ sin 2a (m cos a cos 7 + fc sin a sin 7) , , , , , 

Cg={—, — )^± ^ ^ 3^m,-fc. 2.15 

\dk dm J 2 sin /3 [(fc sin a)2 + (m cos a)2] 2 



One solution of the dispersion relation (2.14) is (w, fc, m) — (1, 1, 1); this is the incident 



wave. In terms of the nondimcnsional variables, the wave fields are 

[ij,b,p] = ^e'(^+^-*)[l,l,-icos/3]+c.c. (2.16) 

where A{t — sz) is the wave envelope and s^^ = — sin2a/2tan/3 is the group velocity in 



the direction normal to the slope (put fc = rn==lin(2.15)) 



2.3. The 'near-criticaV approximation 

We will develop a reductive approximation which is based on taking a distinguished limit 
in which a, /3 — 7, Re~ and Pe~ are all small. To motivate our scaling assumptions. 



we recall the classical solution of the reflection problem given by Phillips (1966). The 



no-flux boundary condition that w = at z = is satisfied by superposing a refiected 



plane wave on top of the incident wave in (2.16). In terms of nondimcnsional variables 



the vertical wavenumber of the refiected wave is 



^^tan(7+^_ 

tan(7 - /3) "■ ' 



The critical reflection condition is that /3 = 7 and then, according to (2.17), the wavenum- 
ber of the reflected wave is infinite. When 7-/3 is small, m is large and the refiected wave 
has a much smaller length scale than that of the incident wave. This singularity, or near 
singularity, signals that the assumptions implicit in the Phillips' solution (stationarity, 
linearity and absence of dissipation) fail when the incidence is near critical. 
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We define the parameter e bym 

, ,1 

e = (a tana) 3 . 



(2.18) 



The reduction is based on the assumption that e is small. As a measure of the degree to 
which the slope departs from the critical condition /? = 7, we introduce a defined 



/3 = 






7 = 






(2.19a, &) 



The sign of a determines if the reflection is supercritical (ct < 0, upslope reflection) or 
subcritic al [a > 0, d owns lope reflection), cr = corresponds to precisely critical reflection. 
From ( ^.17|) and (|2.19|) it follows that when e <C 1, 



tana 



ae^ 



(2.20) 



Thus the 'near-critical' condition in (2.19) ensures that the vertical wavenumber of the 
reflected wave is much greater than that of the incident wave. 

Using (2.15), one can also show that the group velocity of the reflected wave in the 
near-critical case is 

tan 2/3 1 



rhm,l]. 



(2.21) 



tan/3 rn^ 
Using (2.20|), one sees that the z-component of the group velocity in (p. 21) is of order 



£ . Thus if £ ^ 1 it might take an impractically long time to establish the solution given 



by Phillips (1966) 



The theory employs a matched asymptotic expansion in which the incident wave 



in (2.16) is in the outer region and the reflected disturbance is largely confined to an 



inner region, which is essentially a boundary layer close to the slope. 

a useful inner coordinate is likely to be 



From (2.20), we can anticipate that as e 

"^ tana z. The fields of the reflected wave are 



['4',b,p] 



A 

2 



-t) 



sin(/3 -I- 7) 
' sin(/3 - 7) 



icos/3 



(2.22) 



and this suggests the introduction of the more appropriate variables B = te^/sina and 
P = p/cos/3. 

Using ( p. 21 ), we can also estimate the time for the reflected wave to travel through a 
distance on the order of its own wavelength. One finds that this transit time scales as £^ 
which motivates the introduction of a slow timescale t2 ^ e'^t. 

The other scaling assumptions are that 



i?e = 



sm a 



and 



Fe = 



^6 



sm a 

£6 Kg ' 



(2.23a- b) 



The main justification for the choices above is a posteriori — they work in the sense 
that the dissipative terms are comparable to the others in the final amplitude equation. 
Knowledge of s is then the key to recovering a specific physical situation. 

2.4. Nonlinearly reflected second harmonic 



Because of the quadratic terms in ( 2.10 b), one can anticipate, following Thorpe (1987) 
the nonlinear creation of a second harmonic in the small [z = 0(£^)] region in which 
the advective terms become important. However this second harmonic does not remain 

t Because a « cos a, e is bounded as a -^ 7r/2. 
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confined to the region in which z = ©(e^). The second harmonic can radiate into the 
outer region where it appears as a plane wave. A successful completion of the matching 
problem requires that one recognizes this possibility that the inner region drives the outer 
region through this nonlinear process. 

We now assemble the fields which describe the second harmonic in the outer region. A 
modulated plane wave second harmonic with upward group velocity is 



[ip,b,p] = ^H e'^^'=+"''-^''>[l,b,ip] + c.c, 



(2.24) 



where H(t — sz) is the envelope of the second harmonic and s~^ is the vertical group 
velocity. The matching will determine H in terms of the incident amplitude A. In (2.24), 
we use the symbol * to denote quantities associated with the second harmonic. In the 
nearly critical case these quantities are given to leading order by 



. , , 2 a 

m(a) = — — cosec— tana 



COS — + 2 (2 cos a — 1) ' 



6(a) 
p(a) 
S{a) 



1 



m 



COS a, 



(2 cos a — 1)' 



(2.25a) 

(2.25b) 

(2.25c) 

(2.25d) 

If a > 7r/3, the second harmonic is evanescent [j] in z, while for a < 7r/3 the second 
harmonic propagates away from the slope. rh{a) is plotted in figure 0. The angle 9 
between the slope and the direction of propagation of the second harmonic is 



tana (4 + m^ cot^ a) ^ 
2 2 + ifi cot a cot J 



= — arctan 



2 tana 
m(a) 



(a < 7r/3). 



(2.26) 



Both 6 and the group velocity s~^ are shown in figure @. Because both 9 and the group 
velocity are nonzero the second harmonic is not trapped in a boundary layer close to the 
slope. 

Slinn & Kiley (1998a) have reported that when 7 < 7r/6, intrusive layers forms near 
the slope. However, interestingly, in the case 7 = 7r/6 (i.e., the critical value above which 
the nonlinearly reflected second harmonic is evanescent), there is a uniform thickening 
of the dye layer along the slope. 

2.5. The outer region 

The first two ord ers o f the outer solution are obt ained by taking a combination of the 
incident wave in (2.16) and the second harmonic ( p. 24 ) 



['ip,b,p] :^ ^Ae'^''+'-*'^ [l,l,~icosp] + ^H e'^^''+'^'~^'^[l,b,ip] + 0{e'^) + c.c. (2.27) 

The connection between H{t— sz) and A{t — sz) is determined by the matching problem 
in the inner region. 



f Javam et al. (1997) showed the nonlinear creation of evanescent harmonics. Apparently 
they never performed numerical experiments in the case 7 < 7r/6 when they should propagate 
according to the present results. 
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3. The inner region 

In the inner region, we use a stretched coordinate ^ to describe the reflected disturbance 
close to the slope and a slow time scale ^2- Following our earlier discussion in section 2.3, 
these scales are: 

(3.1a, &) 
(3.2a- b) 
(3.3a — c) 



^ = e tana z, ^2 = M^ ^i 

where fi = cot{a/2). As suggested by (2.22), we introduce 



B 



and p = cos [3 P. 



sma 
Finally, it is convenient to define 

u = e^ tana t/, 
so that we have 

U = -^s, and W = *^. (3.4a - h) 

In the zone close to the slope, where <^ — 0(1), B and U are also 0(1). That is, the 
buoyancy perturbation and the along-slope velocity are 0(£^^) larger than in the incident 
wave. 



In the inner region, using ( 2.22 ), (2.9) becomes 



DU 
"Dt 



-P^+B^ 


,£2 


■ a 1 dW 
2^ ^ ot 


U^ + W^-- 


= 0, 





DB 
Ijt 



B = e^ [-fiP^ + iiaB + u^U^^] + 0{e'^), 



U 



where the convective derivative is 

D 
Dt 
Using the complex variable 5* = ^ + iP, (|3.5D can be written compactly as 



(3.5a) 

(3.5b) 
(3.5c) 
(3.5d) 

(3.6) 



S^t - iS'c = -£J{^, S^) +e^ -fJ-S^tz + i^ (f^S^ - S^ 



tJ- 



{Wtt - yWt) 



Y^iBt- iB) + Aig^%c + t^CS^ +Oie% 



where 







Vq +K6 



c 



Uq- Kg 



2/X ' ' 2^ ' 

and J(a, h) is the Jacobian. 

The equation above must be solved with the boundary conditions: 



(3.7) 
(3.8a, 6), 



5(a;,0,t,i2) + 5*(a;,0,t,i2)=0, 
S'4(a;,0,t,i2) + 5'|(a;,0,t,t2)=0, 
%(x, 0, i, t2) - Sl^{x, 0, i, iz) = + 0(e4). 



(no normal flow) (3.9a) 

(no sHp) (3.9b) 

(no flux of buoyancy) (3.9c) 



By taking the inner limit of the right hand side of (2.27), we obtain the matching 
condition 



lim S ^ A{t2 



^— i'OO 



)e'(--') + £i7(f2) tZ£££^e-2'(*-) 
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eH*{h 



, (l+P*sec/3) ^2i(t-x) 



+ 0(e'). 



(3.10) 



The weakly nonlinear analysis proceeds by introducing not only slow space and time 
scales but also seeking a solution of (pj) in terms of the regular perturbation expansions 

S^SQ + eSi+e''S2 + 0{e'^). (3.11) 

Substituting ( |3.11 ) into ( |3.7| ) leads to the following hierarchy: 

e°: (a* -i) 9^^0-0, 

e^: (at~i)9^5i = -J(*o,^oc), 



(3.12a) 
(3.12b) 
(3.12c) 



e^ : {dt - i) d^S2 = - [J(*o, Si^) + J(*i, %)] + [Lr]o. 

In ( 3.12| a), [LTjo means the linear term, in the square bracket of the right hand side of 
(^.7[), evaluated with 5*0. 

3.1. Leading order 
From ( 3.12| a), we obtain the leading order solution 

S-Q = e-'(^^*)5(e, h) + Ae'(^-*) - ^*e-'(^-*\ (3.13) 

and the leading order streamfunction 

*o-^e'(*--)5(^,t2)+c.c.. (3.14) 

The evolution of 5(^,^2) will be determined at higher order. However, at this order, the 
matching condition (3.10) is satisfied provided that 

(3.15) 



lim S 

4— *CXJ 



A* 



whereas the no-flux and no-slip condition condition at ^ = requires that 

5(0, t2) = S^{OM) = %(0,t2) = 0. (3.16) 

Equation ( 3.13| ) is not the most general solution of ( |3.12| a). However, for simplicity we 
harmonic which is required by the matching condition to the incident 

3.2. Order e 



include only e 
wave. 



Equation ( |3.12 b) gives 

Si^t — iS'i^ : 
leading to 

5l=i7^(2;,t,t2)+e2'(*--) 



\ [e2"j(e--5, 6^-54) -f J(e'-5*,e-'-5j)] , 



'-SS( 



S„ du 



S*S^ 



(3.17) 



(3.18) 



We satisfy the no-normal flow condition (3.9a) by requiring that TZ is real. The no-slip 
condition (3.9b) and the no-flux condition (3.9:) are satisfied provided that S satisfies 
these conditions. At this order, the nonlinear effects produce a rectified and a second 
harmonic waves. 

Now we match (3.18) with (3.10). The most important result is that this matching 
condition defines the envelope of the nonlinearly reflected second harmonic: 



/>oo 

Jo 



(3.19) 
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Thus we have an expression for the amplitude of the second harmonic in terms of the 
incident wave. The matching teUs also that TZ = i7e^'^^~*^(l —p sec/3)/2i + c.c. and that 
the rectified flow vanishes in the outer region because of condition (3.15). 



3.3. Orders'^ 
At this order, we have ( 3.12| c). To avoid secular growth, all the resonant terms on the right 
hand side must vanish. This condition determines the evolution equation for 5(^,^2)- 

Although it is not initially obvious, all nonlinear resonant contributions cancel m and 
the final evolution equation is linear 

St,i - iaS^+S- g^S^^^ = A*(t2). (3.20) 

Using ( ^.20 ), one can obtain an alternative expression for the envelope of the second 
harmonic: 



Ht^ + 2iaH = 2q^ 



4d^ 



A\ 



(3.21) 



We emphasize that although there are no nonlinear terms in (3.20), the nonlinearity is 
important for the generation of the second harmonic and also for all the nonlinear con- 
tributions to 5*1 in (3.18). And when we come to visualize the solution, these nonlinearly 
forced components of the solution lead to a symmetry breaking. 

Equation ( ^.20 ) is third order in space and we are imposing four boundary conditions 
in ( 3.15D and (pT6 ) . Consequently, the problem is overspecified and we resolve this issue 
by discarding the no buoyancy flux condition 5^^ = 0. This unsatisfactory point might 
be corrected by demanding that the Pcclct number in (2.11b) be very large. In this 
circumstance, one would expect a very thin buoyancy diffusive layer in which dynamics 



similar to that of Wunsch (1970) and Phillips (1970) is important 



4. The inviscid case 



We first c onsid er the special case of ( 3.20 ) in which the fluid is inviscid {g — 0). The 
solution of (3.20) which satisfies the initial condition that 5(^, 0) = is 



S{iM) = f A*{t2 - r)e'^^\ll Ji (2Ve^) dr, 



(4.1) 



and (3.21) leads to the following expression for the envelope of the second harmonic: 



H{t2 



g2i.(r-t.) ^2(^) j^_ 



(4.2) 



The amplitude ^(^2) of the incident wave must be specified to completely determine 
the solution. As in all the experiments, we take the simplest case in which A{t2) switches 
on suddenlyra at ^2 = 0. That is, A{t2) = 1 when ^2 > 0. 



4.1. The critical case 



In the critical case a — Q the integral in (4.1) can be evaluated and 

5 = 1- Jo(2v/et^), 



(4.3) 



I The underlying reason for the miraculous cancellation of the resonant nonlinear terms is 
the following special case of the Jacobi identity; 

J [Q, J{Qi, Q*)l + J [Q\J{Q, Qe)] + J [J[Q, Q*), Qe] = 0. 

:|: Because A depends only on the slow time i2, the sudden switch-on means that the incident 
wave achieves its ultimate constant amplitude on a time scale which is slow relative to t but 
fast relative to i2. 
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where Jo is the Bessel hmction of the first kind of order 0. Thus, in this critical case, the 
solution has a 'similarity' form in which the thickness of the inner region is inversely pro- 
portional to time. That is, there is no steady solution as t — > oo. Instead, the disturbance 
near the slope becomes strongly oscillatory as the undulations of Jq are intensified. 



Using (4.3) we have: 
[*o,C^o] =cos(x-t) 



[Wo,Po,Bo]^siii{x-t) 
'Ji 



l-Jo,-2h — 
X. 



Jo-l,l + Jo,-2t2 — 



*1 

Ui 



t2 



2ti 



X 



(1 - Jo) + Jo^ + Ji - 1 cos2(x - t) + t2(l - Jo) 



J2(l - Jo 



Wi - 2i2 1 



X 

j2 _ j2 



j2\ /J 

-^\cos2{x-t) + 2tl f ^(1- Jo) 

Jl 



•II 

X 

^2 



X 



(1- Jo) sin2(x-i). 



Bi — 2^2 



4(l-Jo) + 4)sin2(x-t), 



X 



X' 



(4.4a) 
(4.4b) 
(4.4c) 
(4.4d) 
(4.4e) 
(4.4f) 



where x = 2\/p2 and Ji = Ji (x) . Notice that in this critical case the along-slope velocity 
Uq and the buoyancy perturbation Bq both grow linearly with time. This response is 
analogous to that of a resonantly forced oscillator. 

The streamfunction is shown in figure g. At small times the reflection process creates 
a regular array of counter-rotating vortices. As time progresses, figure ||(b) shows that 
the scale of the vortices decreases. Panels (c) and (d) are both at TV sin /3 t = 3. The 
"velocity vector" presentation in panel (d) more clearly displays the asymmetry of the 
vortices which is the effect of the nonlinear terms in \E'i. 

Figure ^ shows the distortion of the isopycnals as the oscillations amplify. In panel (a) 
the disturbance is very small and one sees essentially the initial background stratification. 
In panel (b) the disturbance begins to 'fold-up' the isopycnals and, near the slope, this 
process produces a region of static instability. In the panel (c) the development of small 
scales in the isopycnal field is evident. 

In this case, with cr = 0, the amplitude of the second harmonic is a linear function of 
time: 

^2 



H{t2) ^-h = - 



tan 7 



t, 



(4.5) 



and so e^i becomes comparable to ^o when ^2 — 0(e^^). For these reasons the expansion 
becomes disordered when ^2 = ^'(e^^) and the results above are no longer reliable. 
However, as we show below, well before this breakdown, the buoyancy becomes statically 
unstable. Thus the expansion above strongly suggests that the next evolutionary stage 
is characterized by the onset of turbulence triggered by overturning instability. 

4.2. Overturning instability or stratified shear flow instability? 

One scenario for the transition to turbulence is that the growing disturbance produces 
a statically unstable density field which then overturns (we refer this as overturning in- 
stability) . An alternative is that the local Richardson number might fall below a critical 
value while the density is still statically stable. In this second case the overturns are 



produced by a rapidly growing secondary shear flow instability. Experimentally, Thorpe 



fe Haines (1987) have reported that the overturning instability is very likely to be con- 
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vectively driven. On the other hand, Shnn & Riley (1998a) identified a shear instabihty 



mechanism; but they noted that the Reynolds number of their numerical simulations was 
matched to experimental values by forcing larger amplitude waves in a more viscous fluid. 
Consequently, the Richardson number in these simulations was relatively low compared 
to those of laboratory experiments. However, all these claims must be viewed cautiously 
because both instabilities are intrinsically related in a stratified flow. 

Using our analytic results we can make a rough assessment of these two possibili- 
ties. Taking into account the background linear buoyancy, and the first order correction, 
we can calculate the overturning time to which is the time at which negative vertical 
buoyancy gradients first occur. In dimensional variables: 

-v/8tan7 3 , , „, 

to= l, . / e-^. 4.6 

N sm 27 

We can also calculate the time tg at which the minimum Richardson number first falls 
below 1/4: 



ts^\ to>to. (4.7) 

V cos 7 

This analytical calculation also shows that both of the unstable conditions above occur 

first at the wall (see also Javam et al. (1997)| ). 



Because to < tg we can argue that the convectively driven overturning instability 
should appear first. In the typical experimental case e — 0.3, for example, both critical 
times are plotted as a function of the slope angle 7 in figure 0. It is interesting to note 
that N suvfS to is an increasing function of the slope angle 7. 

However, let us note that the domain of validity of the Miles-Howard theorem (sta- 
bility if Ri > 1/4) applies to steady, parallel shear flows, whereas the present flow is 



unsteady and nonparallel. In addition, Eq. (4.7) was derived using the usual definition of 
the Richardson number Ri — N'^ /U"^ (see, for example, Kundu (1990)[ ). An alternative 
definition could be 

ffi=-^^^, (4.8) 

where the numerator takes into account not only the background buoyancy, but also the 



perturbative part of the buoyancy. It is clear from Eq. (4.8), that Ri would reach 1/4 
before zero, and therefore, before negative vertical buoyancy gradients first occur. So, 
the whole discussion simply indicate that there are plausible reasons to expect that the 
solution we have found will become unstable, and ultimately turbulent, and we have three 
criteria for instability; the static instability time is sandwiched between the Ri — 1/4 and 
Ri = 1/4 times. 

4.3. The noncritical case 

We now turn to the case in which the incident wave is not precisely critical so that a 7^ 0. 
We continue with the assumption that the fiuid is inviscid {g = 0) and that the incident 
wave envelope is ^(^2) = 1 if ^2 > 0. 

In this case (3.20|) has a steady solution which satisfies the boundary condition at 



? = 0: 

5=1- e-'i/". (4.9) 

This steady solution is an approximate version of the linearly reflected wave identifled 



by Phillips (1966). Notice how the scale of the oscillations is reduced as ct ^ 0. Thus, 



when the incidence is nearly critical, one expects to see an initial reduction in scale 
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which is the t ^ behaviour identified in (3.2C). But this scale reduction is arrested at 



time t2 — 0{l/\a\) when the Philhps solution in (4.9) is established as a steady state. 

In order to understand the details of how the steady solution in (4.9) emerges we can 
use the solution of the initial value problem given in (4.1): 



5 = / e''^Ji{u) du, 

Jo 



[*o,Wo] = [Re,Im] 



2V?*2 



Mu) e'(*+^-") dt 



[Uo^Bo 



-Re, Im] 



JoiW^r 



Ut + CTT — X) 



dr, 



(4.10a) 

(4.10b) 
(4.10c) 



[Ui, Bi] = -[-Rc,Im] 



^2i{t-x) 



{SS^^ - si) + {S*Si)^ 



H 



-a sin crio 



(4.10d) 
(4.10e) 

At early times, the streamfunction in figure ^ is similar to the critical case in figure ^. 
However, as time progresses, the steady solution in (19) is set up first in the neighbor- 
hood of the wall, and then this cellular pattern expands outwards. Figure ^ shows the 
region in which the flow becomes steady is characterized by regular pattern of vortices. 
Figure Uu shows the isopycnals at a given time: the pattern is similar to experimental 
results reported in figure 12 of Thorpe (1987) (see also Mac Phee 1998). 

As in the critical case, one can compute to and tg for different values of a. Again, the 
instability is convectively driven and initiated at the wall, to is plotted in figure O as a 
function of the slope angle 7 for different values of a. An interesting point is that, for a 
given value of 7, the positive values of a leads to an earlier appearance of the instability. 
This point is rationalized by noticing th at, o nce 7 is given, a positive a corresponds to 
a smaller value for a. As the expression (4.6) is an increasing value of the angle a = 27, 
near-critical up-slope (respectively down-slope) reflections are unstable slightly before 
(respectively after) critical reflections. Indeed, plotted as a function of a, to and ts are 
almost independent of a. 

that for 



This is consistent with the observations reported by De Silva et al. (1997) 



moderately supercritical waves the instabilities developed near the bed. They have, how- 
ever, also studied experimentally the variation of the boundary layer thickness as the 
incident waves become far from critical, and their results show that the instability is ini- 
tiated away from the bed. In the framework of this near-critical reflection theory (i.e. for 
small values of e and a), we found that the wave overturning always starts on the slope, 
however as time continues the unstable region extends away from the boundary (for ex- 
ample, see the unstable region in flgure 10). For strongly subcritical and supercritical 



cases, internal wave reflection from the sloping bed should be interpreted as w ave- wave 
interaction b etween the incident and the reflected waves since, as shown by De Silva 
et al. (1997) , the area of interaction region increases progressively as the waves depart 



from critical condition. 



5. Viscous effects 

In the viscous case, g ^ 0, there is a steady solution of ( 3.20| ) even if the forcing is 
precisely critical. For the sake of simplicity, consider the case in which g y^ 0, a = 0, and 
the incident wave envelope is ^(^2) = 1 if ^2 > (this is a typical experimental switch 
on). The steady solution of ( |3.20 ) which satisfies the no- mass flux and no-slip boundary 
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conditions at ^ = is analogous to the western meridional boundary layer, also called 
the Munk layer (see for example Pedlosky 1987). It reads 



(5.1) 



(we continue to assume that k ~ and discard the no buoyancy flux boundary condition 
at ^ = 0). The solution above, presented in figure |l2| in the extremely viscous case i = 1, 
shows that 20^'^ is the viscous boundary layer thickness in the terms of ^; in dimensional 
variables this boundary layer thickness is 

1/3 / „ \ 1/3 

(5.2) 



V + K 



1 



COS a 



Thus, in the steady state, the viscous boundary layer thickness is a decreasing function 
of the slope angle. 

The steady state solution is probably irrelevant in many experimental systems because 
g ^ 1 (see table |^) and one expects that turbulent transition occurs before the steady 
state is approached. In this case with p <C 1, we can again use asymptotic matching to 
develop an approximate solution of the initial value problem. There is an interior region 
in which the effect of viscosity are small and the solution is approximated by m^- 
However this interior solution of section H does not satisfy the no-slip condition and so 
it is necessary to include a viscous sublayer close to the topography (see figure ^ for a 
schematic representation of the different regions) . 

The details of this matching problem are in the appendix. The solution in the interior 
region is 



S{tt2) = 



t2 



A*(i2-r)e" 

t2 



I Ji ( 2^^T ) dr 



-g I {S^,it2-T)-ic7Su{t2-T))Jo{2^^T) dT + 0{g') (5.3) 



Si4t2) = -^ / Vh^^ e'^(*^-^) A*{T)dT 



where 





and Jo is the Bessel function of the first kind of order 0. 
In the viscous sublayer the solution is 

5(e,t2)-e /'e'"'*^""^A*(T)dT 

Jo 



(5.4) 



+Q 



SUtt2)+ I du- 





-i^l^e^u i't2~y. 



A*(r)e''"(*^"^M7 



+ 0{g'). (5.5) 



In the critical case with a = and A{t2) = 1 if t2 > 1 the integrals above can be 
simplified. The interior solution is 



S{^,t2)^l-Joi2V^h) 



2^^ 
and in the viscous sublayer, the solution is 



sin2V^ - 2V^ cos 2y^W 0(^2) (5.6) 



S{iM)=£.t2- 



gt 



du- 



-«V4e^t2« 



Oig"). 



(5.7) 
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Both functions are plotted against z in the figure n^ in the typical case p = 0.1 and 
e = 0.3. It is clear that the region where the viscous effects are important corresponds 
only to a very thin region along the slope. 

6. Conclusion and discussion 

The thrust of this paper has been to study the weakly nonlinear and nearly critical 
incidence of internal waves onto a slope. The scalings of sections 2 and 3 amount to taking 
the distinguished limit |/3 — 7I -^ with a ex |/3 — 7^/^ (a is the amplitude parameter 
in (2.12)). At leading order, these assumptions give the linear oscillator equation (3.12a) 
in which the coordinate normal to the slope, ^, appears only parametrically. Thus, the 
buoyancy oscillations along the slope are uncoupled at leading order. The weak coupling 
between oscillations at different ^'s is uncovered by higher orders in the expansion scheme 
and is apparent in the forced dispersive wave equation (3.20). One can then view the 
incident internal wave as a nearly resonant forcing of this continuum of weakly coupled 
(and weakly damped if p ^ 0) oscillators. 

The scenario above describes the initial evolutionary stages of nearly critical inci- 
dence. However the limitations of this approach become apparent when the oscillations 
become so extreme as to either overturn the buoyancy field, or strongly violate the Miles- 
Howard stability condition (see figure 6(c)). In either case, we expect a rapid transition 
to turbulence, dramatically enhanced mixing in the neighbourhood of the slope, and the 
production of intrusive layers (e.g., Ivey & Nokes 1989; De Silva 1997; Mac Phee 1998). 

Within the present framework the most interesting complication which can be included 
is oblique incidence. Experimental data by Eriksen (1998) and theoretical work by [Thorpe 



(1997) have recently emphasized the importance of alongslope currents in the refiection 
process for obliquely incident waves in a uniformly stratified rotating fluid. We speculate 
that the weakly nonlinear term will have interesting consequences, such as mean flow 
induction, if the incidence is oblique. 
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Appendix A. Derivation of the solution in the viscous case. 

Defining 

J{S,M)^S{S,M)e~'"'- and g{h) ^ A* {t^) e'"''-, {Ala,b) 

( [3.20 ) becomes 

/t.? + /-e'/«C-5 (A 2) 



with the 4 conditions 

/(0,i2) = /c(0,t2) = /«(0,i2) = and lim 7(^,^2) = sfe). {Aia-d) 

4 — >oo 

A.l. Interior region 



In the interior region, the solution of (A 2) is obtained expanding / in powers of g: 

fiiM) = MiM + Q h{£.M)+0{g^). (A 4) 
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The substitution of (A 4) into (A 2) leads to the foUowing hierarchy: 

0° ■ /otaC + /o = g, 

g' ■■ /it.c + /i - 0. 

At leading order, p°, the solution of (A 5a) should satisfy the 2 conditions: 



/o(0,t2)-0 



and 



lim fo{Lt2)^g{t2). 



9it2^r) J^J,{2^^r dr 



The solution can be obtained using the Laplace transform and we get 

Jo 
where Ji is the Bessel function of the first kind of order 1 
At order g, (A 5b) with the conditions 

/i(0,i2) = /i.(i2) and 
leads to the solution 



/i(e,t2)=/ /i.(i2 - r) Jo (2v/^) 



hm /i(e,i2) = 0, 



dr. 



(A 5a) 
(A 5b) 

iA6a,b) 

(A 7) 

{A8a,b) 
(A 9) 



The unknown function /i*(i2) will be determined by matching the viscous sublayer so- 
lution. 

A. 2. Viscous sublayer 

I n th e viscous sublayer, using the stretched coordinate rj — ^/ g and h{ri, ^2) = f{(,, t2)l Q, 
(A 2) becomes 

ht2ri - hr,,,,, + g h = g. (A 10) 

Expanding h in powers of g, the leading order solution is 



hoiv, h) = fi*{t2) +V 9{T)dT + r(?7, t2), 

JQ 



where r (77,^2) satisfies the diffusion equation 
with the three conditions 



(All) 
(A 12) 



r^{0,h) 



5(r)dT, lim r = and ?'(0, ^2) = ^ /i*(^2)- {A13a — c) 

^ — >CXD 



Let us stress that the two first conditions define completely the solution of ( A 12| ), whereas 
the last one determines /i*(t2) in terms of g{t2)- 

Using the Laplace transform, one obtains the solution of ( A 12| ) and we finally get, not 
only the function 



flAt2 



'T^ Jo 



^/h - T 5(r)dT 



but also the the general solution in the viscous sublayer given by 



fiU2) = U 9{r)dT + g 
Jo 



h*{t2)+ I du 





t2 g-C /4ru /•t2~U 



g{T)dT 



(A 14) 



+ 0{g'). (A 15) 
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Figure 1. Schematic view of the reflection of the incident wave. The angle between the bottom 
slope and the horizontal is 7; the angle between the incident group velocity and the horizontal 
is l3, and a = /3 + 7. Cj indicates the group velocity and g indicates gravity. 
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Second 
Harmonic 




Figure 2. Schematic view of the reflection process for 7 ~ /3. The reflected wave is alongslope, 
whereas the nonlinearly reflected second harmonic, represented with the dotted line, makes an 
angle 6 with the slope. 
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Figure 3. s^^ (solid line), the 2:-component of the group velocity, and Q (dashed line) angle 
of reflection between the slope and the direction of propagation of the second harmonic in the 
critical case. 
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Figure 4. m(a), the z-wavenumber of the nonhnearly reflected second harmonic, in the 
near-critical case. When a > n/3, m{a) is complex. 
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Figure 5. Streamfunction for four different times in tfie critical case. A^ sin/3 t — ^, 1 and 3 in 
the panels (a), (b) and (c). Panel (d) is an alternative "velocity vector" visualization of panel (c). 
e = 0.3, y = f} — 20°. In panel (d), there are two pairs of counter-rotating vortices immediately 
adjacent to the slope. The clockwise vortices are slower and thinner than the counter-clockwise 
vortices. This symmetry breaking is a result of the second harmonic term such as >l/i in (4.4). 
The dimensions of the panel are 2A/sina in the a;-direction (i.e. two alongslope wavelengths) 
by A/sinQ in the z-direction. 
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Figure 6. Buoyancy field in tiie critical case at Nsin/S t —3, 10 and 20 in the panels (a), (b) 
and (c) respectively, e = 0.3; /3 = 7 = 20°. The dimensions of the panel are A/ sin a in the 
a:-direction by A/5 sin a in the ^-direction. In the panels (b) and (c), some overturned regions 
are evident. 
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Figure 7. NsinP to (dashed line) and N sin f3 ts (solid line) as a function of the slope angle 7. 
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Figure 8. Streamfunction for four different times in the near-critical case: /3 = 22.58°; 
7 = 17.42°; a = ~1; £ = 0.3. iV sin/3 t = 5, 10, 15 and 20 in the panels (a), (b), (c) and 
(d). The dimensions of the panel are A/ sin a in the a;-direction by A/2 sin a in the z-direction. 
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Figure 9. Streamfunction in the near-critical case at TV sin/? i = 100: /3 = 22.58°; 7 = 17.42°; 
a = —1; e — 0.3. The dimensions of the panel are A/ sin a in the x-direction by A/2 sin a in the 
2-direction. 
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Figure 10. Buoyancy field at Nsm/S t = 20 in the near critical case /3 = 22.58°; 7 = 17.42° 
and e — 0.3. The dimensions of the panel are A/sina in the a;-direction by A/5 sin q in the 
2-direction. 
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Figure 11. N &vaj5 to against the slope angle 7. a = (solid curve), (7 = 1 (dotted curve), 

a = — 1 (dashed curve). 




Figure 12. Final steady state solution in the critical case: /? = 32.58 ; 7 = 27.42 ; a = —1; 
e = 0.3. The dimensions of the panel are A/ sin a in the x-direction by A/4 sin a in the a-direction. 
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Figure 13. Schematic view of the different regions used in the matched asymptotic expansion. 
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Figure 14. Function / against the distance normal to the slope z at time t — 1 and for q = 0.1, 
in the critical case. In the inset, we plot the viscous sublayer solution valid only very close to 
the wall, i.e. 2 <C 1. 



